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ABSTRACT 

The discovery of many planets using the Kepler telescope includes ten planets orbiting eight binary 
stars. Three binaries, Kepler-16, Kepler-47, and Kepler-453, have at least one planet in the cir- 
cumbinary habitable-zone (BHZ). We constrain the level of high-energy radiation and the plasma 
environment in the BHZ of these systems. With this aim, BHZ limits in these Kepler binaries are 
calculated as a function of time, and the habitability lifetimes are estimated for hypothetical ter¬ 
restrial planets and/or moons within the BHZ. With the time-dependent BHZ limits established, a 
self-consistent model is developed describing the evolution of stellar activity and radiation proper¬ 
ties as proxies for stellar aggression towards planetary atmospheres. Modeling binary stellar rotation 
evolution, including the effect of tidal interaction between stars in binaries is key to establishing the 
environment around these systems. We find that Kepler-16 and its binary analogs provide a plasma 
environment favorable for the survival of atmospheres of putative Mars-sized planets and exomoons. 

Tides have modified the rotation of the stars in Kepler-47 making its radiation environment less 
harsh in comparison to the solar system. This is a good example of the mechanism first proposed 
by Mason et al. Kepler-453 has an environment similar to that of the solar system with slightly 
better than Earth radiation conditions at the inner edge of the BHZ. These results can be reproduced 
and even reparametrized as stellar evolution and binary tidal models progress, using our online tool 
http://bhmcalc.net. 

Keywords: binaries: general - planet-star interactions - planets and satellites: individual (Kepler-16b, 
Kepler-47c, Kepler-453b) - stars: activity 


1. INTRODUCTION 

The outstanding work of the Kepler Observatory cir- 
cumbinary team has resulted in the discovery of 10 con¬ 
firmed planets around eight moderately separated bina¬ 
ries (all have binary periods in the range of 8-60 days). 
These are Kepler-16 (Doyle et al. 2011), Kepler-34 and 
Kepler-35 (Welsh et al. 2012), Kepler-38 (Orosz et al. 
2012b), Kepler-47 (Orosz et al. 2012a), Kepler-64 (Kos- 
tov et al. 2013; Schwamb et al. 2013), Kepler-413 (Kostov 
et al. 2014), Kepler-47 (?) and Kepler-453 (Welsh et al. 
2015). Interesting, three out of eight of these systems 
(38%), namely, Kepler-16, Kepler-47, and Kepler-453, 
host giant planets in the circumbinary habitable-zone 
(BHZ)% 

The study of conditions that planets face while orbiting 
a binary system has attracted the attention of the exo¬ 
planetary community. Several authors have studied the 
stability of their orbits (Dvorak 1986; Holman & Wiegert 
1999) and the kind and level of insolation experienced by 
planets while illuminated by two stars (Harrington 1977; 
Haghighipour & Kaltenegger 2013). The conditions for 
the formation and migration of these planets have been 
studied by Pierens & Nelson (2008), Kley & Haghigh¬ 
ipour (2014), and Kley & Haghighipour (2015). Efforts 
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are being made to constrain the magnetic and plasma en¬ 
vironment of circumbinary planets (Mason & Clark 2012; 
Mason et al. 2013; Sanz-Eorcada et al. 2014), finding that 
binaries could potentially pose severe restrictions on the 
habitability of Earth-like planets and exomoons. We now 
know that stellar aggression, in the form of high-energy 
radiation and coronal-activity-related mass-loss, would 
have an effect on the evolution of planetary atmospheres, 
especially their capacity to retain water (see, e.g. Lam- 
mer 2013). 

In Mason et al. (2013) we proposed a mechanism for 
which some BHZ planets experience Earth-like or lower 
levels of stellar aggression and thereby have the potential 
to retain water. Tidal torquing of stellar rotation, espe¬ 
cially for certain binary periods and initial eccentricities, 
aids in the reduction of stellar activity. The reduction 
in intensity of rotationally driven dynamo action results 
in a decrease in stellar coronal activity, potentially pro¬ 
moting life on planets. We refer to the sum of these and 
other beneficial effects for life on planets around binaries 
as the Binary Habitability Meehanism (BHM). 

Recent observations lend support to the tidal astro¬ 
physics of BHM. Eor example, Wasp-18, a 0.63 Gyr old 
E6 star, is orbited by a nearby giant planet, with a 0.94 
day orbital period. The star shows a low level of activ¬ 
ity, normally expected for a much older star (Pillitteri 
et al. 2014). This effect has been described as “prema¬ 
ture aging.” In another example, the interaction of HD 
189733 with a close-in planet, has resulted in the oppo¬ 
site effect on a much larger timescale (WoIk et al. 2011). 
HD 189733, a relatively old star, displays anomalously 
high activity levels, lending support to our prediction 
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that early tidal rotational torquing could result in a high 
level of activity. We call the latter phenomenon the “for¬ 
ever young” effect, in reference to the fact that tidal in¬ 
teraction of a star with a stellar or substellar companion 
could make the star very active (to appear young) for 
a longer time than expected from single-star evolution 
(Mason et ah 2013; Sanz-Forcada et ah 2014). 

Although the specific physical mechanisms operating 
in the case of Wasp-18 and HD 189733 may be differ¬ 
ent from those predicted in stellar binary systems, these 
discoveries confirm the necessity of including tidal inter¬ 
action in order to assess stellar activity levels of stars 
with companions and its effect on their planets. 

In previous work, we used first order estimation of tidal 
synchronization times (Mason et al. 2013). More recently 
(Mason et al. 2015), we explored habitable niches around 
hypothetical circumbinary environments, using a basic 
model for stellar rotation evolution including stellar tides 
but applied only to the early phases of main-sequence 
evolution. 

Here our models are refined by (1) extending the time 
frame to the critical pre-main-sequence phase where ro¬ 
tational and activity evolution is much richer; (2) ap¬ 
plying a more physically motivated model for stellar ro¬ 
tational evolution, including but not restricted to angu¬ 
lar momentum (AM) transport inside the star; and (3) 
estimating mass loss and X-ray emission using the lat¬ 
est solar-inspired models, relating basic stellar proper¬ 
ties and rotation to chromospheric and coronal activity 
(Cranmer et al. 2013). 

For these three critical refinements, we are applying 
models that have been developed and tested in the case of 
single-stars. We are assuming here that the same models 
can be applied to describe low-mass stars in moderately 
separated binaries (orbital periods larger than 8 days and 
separations larger than several tens of stellar radii). We 
will justify and test this assumption later on in the paper. 

We apply our improved and extended model to the 
interesting and well-characterized Kepler binaries with 
BHZ planets. The planets currently known in those sys¬ 
tems are Neptune to Saturn sized, so no surface habit¬ 
ability is expected. However, it is reasonable to suspect 
that Earth-like planets may exist in the BHZ of similar 
binaries (Eggl et al. 2013). Additionally, since circumbi¬ 
nary giant planets appear to be common, some of them 
may harbor potentially habitable exomoons (Heller et al. 
2014). These planets and putative moons are expected 
to experience harsh environments from aggressive host 
stars. Thus, understanding the radiation and plasma 
environment experienced by these bodies is a key effort 
to accessing their potential habitability (Zuluaga et al. 
2012; Heller & Zuluap 2013). 

The present aim is to constrain the radiation and 
plasma environment of these HZs and to investigate the 
effect that these conditions have on the atmospheric evo¬ 
lution of hypothetical terrestrial planets and exomoons. 

This model can readily be applied to other Kepler bi¬ 
naries or any other system with well-determined param¬ 
eters. Euture telescopes like PLATO (Rauer et al. 2013) 
will likely discover circumbinary planets exhibiting some 
potential for habitability. Our aim here is to inform such 
searches by showing how a more comprehensive theory 
of planetary habitability can be applied to circumbinary 
planets especially in the presence of interacting stars. 


In order to provide the reader an opportunity to repro¬ 
duce these results and to further explore the vast param¬ 
eter space of circumbinary planets, a Web-based BHM 
Caleulator^ BHMcalc , is made available using the link 
http://bhmcalc.net. In Mason et al. (2015) we intro¬ 
duced the first version of the tool. Here we present a 
comprehensive version of the calculator, detailing all the 
new physical effects described in this paper, and provid¬ 
ing an improved and user-friendly interface. 

Other online tools have recently been made available 
(Cuntz & Bruntz 2015; Muller & Haghighipour 2014). 
Like these, the BHMcalc provides renditions of the in¬ 
stantaneous and continuous BHZ. In addition, BHMcalc 
calculates time-dependent stellar and planetary proper¬ 
ties and utilizes them to study HZ environments. Most 
critically, the rotational evolution of the stellar compo¬ 
nents is derived, allowing for estimates of the resulting 
stellar mass loss and magnetic activity. Stellar emis¬ 
sion is used to evaluate effects experienced by planets. 
Specifically, we select the integrated X-ray and extreme- 
ultraviolet radiation (XUV) and stellar wind (SW) fluxes 
as derived for planets, located over a range of distances 
from the binary center of mass, as proxies for habitability. 
Eor this purpose we estimate the XUV- and SW-induced 
planetary atmospheric mass loss. In all cases, the evolu¬ 
tion of each star is calculated independently and orbital 
averages are used to calculate the combined effect of the 
binary on the planet. 

This paper is organized as follows: In Section 2, the 
sample of three Kepler binaries possessing BHZ planets 
is described. Section 3 presents an efficient model to cal¬ 
culate BHZs and their evolution in time. In Section 4, we 
revisit the argument for the need to access the radiation 
and plasma environment of planets as potential drivers of 
atmospheric evolution. High-energy radiation and winds 
are key to constraining habitability. The comprehensive 
model used to calculate the evolution of rotation and 
activity is described in detail in Section 5. The model 
is applied to constrain the radiation and plasma envi¬ 
ronments of the solar system as a reference in Section 
6. The binary rotational evolution affects potential cir¬ 
cumbinary worlds and so is presented in Section 7, and 
the radiation and plasma environments are investigated 
in Section 8. Einally, in Section 9, a summary and the 
conclusions of this investigation are presented. 

2. KEPLER CIRCUMBINARY HABITABLE-ZONE PLANETS 

It is interesting to note that while none of the currently 
known transiting circumbinary planets are Earth-like or 
even super-Earths, this may be due purely to three se¬ 
lection effects. The Kepler telescope is most sensitive to 
the detection of large planets, orbiting small stars, on 
short-period orbits. However, the fact that on average 
circumbinary planets have longer periods than those Ke¬ 
pler planets with single-star hosts is not a selection effect. 
Short-period circumbinary planets are limited by dynam¬ 
ical stability (Harrington 1977; Dvorak 1986; Holman & 
Wiegert 1999) and probably also because of formation 
(Pierens & Nelson 2008; Kley & Haghighipour 2014). 
Remarkably, the BHZ extends well beyond the orbital 
stability limit in many cases (see Mason et al. 2015). In 
addition, residence of planets in the BHZ appears to be 
common, as 3 of the 10 known transiting planets are lo¬ 
cated in the BHZ and three binaries out of eight harbor 
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a BHZ planet. 

In Table 1, properties of the three transiting BHZ plan¬ 
ets and their host stars are listed. Hereafter, we will 
refer to this as the KBHZ sample. The binaries in the 
KBHZ sample are among the best-characterized binaries 
known, mainly owing to the extraordinary observational 
constraints provided by transiting circumbinary planets. 

To study the binaries in the KBHZ sample^ we focus 
on three different scenarios, selected to answer the fol¬ 
lowing questions: (1) If the actual circumbinary plan¬ 
ets were Earth-like planets instead of giant planets, then 
would those planets be habitable?^. From necessity, this 
must be based on current understanding of Earth-like 
habitability. (2) Assuming that the actual circumbinary 
planets, which have masses larger than or around that of 
Saturn, would be able to harbor a Mars-sized exomoon, 
then would these moons be able to retain a dense and 
moist atmosphere in order to sustain life? And maybe 
most importantly, binaries like these are very common, so 
(3) what are the constraints on habitability of currently 
unseen planets within the BHZ of similar binaries? In¬ 
sight into these and other questions is gained by studying 
well-observed cases, where the environmental conditions 
experienced by planets under the reign of the binary may 
be investigated. 

It is important to stress that the current investiga¬ 
tion is not restricted by dynamical or formation con¬ 
straints, which may or may not limit the likelihood of 
these scenarios (Chavez et al. 2013, 2015; Andrade-Ines 
& Michtchenko 2014; Forgan et al. 2015; Pilat-Lohinger 
et al. 2014). The basic dynamical constraint imposed 
planets have stable orbits is included as it extends into 
the BHZ in some cases (see Section 3). We assume that 
exomoons and multiple planets could lie in stable orbits 
within the BHZ of these binaries, based on our own basic 
tests using numerical orbital integrations. 

3. THE CIRCUMBINARY HABITABLE-ZONE 

Estimates of BHZ boundaries quickly followed the Ae- 
p/er discoveries (Mason & Clark 2012; Quarles et al. 2012; 
Haghighipour & Kaltenegger 2013; Kane & Hinkel 2013; 
Mason et al. 2013). Our approach, as that of others, is to 
use the models of Kasting et al. (1993), updated by Kop- 
parapu et al. (2013b) derived for single stars to also esti¬ 
mate the more complex BHZ limits. An improved version 
of the first-order method presented in Mason et al. (2013) 
is developed and applied here. 

HZ limits are determined by the critical fluxes at which 
planetary surface temperatures are compatible with the 
existence of liquid water. Using one-dimensional atmo¬ 
spheric models, Kopparapuet al. (2013b) calculated crit¬ 
ical fluxes assuming that the planet is exposed to black- 
body radiation with an effective temperature Tgff: 

Sef^,i = S'eff,!,© + UiT* + hiT^ “h c{T^ + d\T^ (1) 

Here T* = Teff — 5780A, and 5'eff,i and ai,6i,Ci,di are 
the effective critical flux for a solar-mass star and the 
interpolation coefficients for the ith boundary (see Table 
3 in Kopparapu et al. 2013a). 

^ It is important to stress here that this is an hypothetical sce¬ 
nario. We are not under any circumstance assuming that the dis¬ 
covered giants planets can be habitable or that terrestrial planets 
also exist in the BHZ of all of them. 


The limits of the HZ, around either single-stars or a 
binary, are calculated by finding the radius of a circle 
centered in the barycenter of the system, where the av¬ 
erage flux {S) received from the star or the binary is 
equal to the critical flux calculated with Eq. (1). 

In the case of a single-star, {S) is equal to the nearly 
constant flux iSsing = LjcP measured at distance d. Here 
L is the bolometric luminosity of the star measured in 
solar units and d is in astronomical units (AU). More 
precisely and as we will explain later, L is a corrected 
luminosity that takes into account the response of the 
atmosphere to the incident stellar radiation (see Eq. 4). 

The radius of the Ah HZ edge /ging i will then be given 
by: 

5'efr,i = J2 (2) 

^sing,i 

Around binaries, planetary atmospheres, even if plan¬ 
ets are in a circular orbit, will be subject to a variable 
flux and spectrum of radiation. The instantaneous bolo¬ 
metric flux on a circumbinary planet is calculated as 

" Riile) + Riile) 

where Ri and i?2 are the distances from the stars to 
a point in a circle of radius d. 0 is the angle between 
the line joining the stars and a fixed point on the circle. 
This angle is a function of time, and it will depend on the 
relative angular velocity of the binary and a test particle 
on the circle. 

Ri and R 2 are given by 

RI {d,0) = {d — r 2 sin 6>)^ cos^ 0 

Rl{d, 6>) = (d + ri sind)^ -h rl cos^ 0 

where ri = qr/{l-\-q) and r 2 = qvi are the instantaneous 
distances of the stars to the barycenter and q = M 2 /Mi 
is the binary mass ratio. The relative distance r between 
the stars is also a function of time for the general case of 
an eccentric binary orbit. 

Since the two stars have different effective tempera¬ 
tures, 5'bin and cannot be directly compared as in 
the case of single-stars. However, in the case of stellar 
twins, the two sources will have the same effective tem¬ 
perature and their fluxes can be summed directly. On 
the other hand, in disparate binaries the ratio of the flux 
of the secondary to that of the primary scales with the 
fourth power of the ratio of their effective temperatures. 
As a result the combined spectrum, to first order, may be 
assumed to have an effective temperature equal to that of 
the primary and the total flux equals 5'bin (Mason et al. 
2013). A verification of this approximation is discussed 
below. 

A more consistent method to take into account the 
differences in effective temperatures of the stellar com¬ 
ponents is to replace bolometric luminosities Li and L 2 
in Eq. (3) by weighted luminosities (Haghighipour & 
Kaltenegger 2013; Cuntz & Bruntz 2015): 

L© = Li, 2(1 + [-5eff.i(Ti,2) - (4) 

The weights in parenthesis account for the different 
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Parameter 

Kepler-16b 

Kepler-47c 

Kepler-453b 

Stellar Properties 

Ml (Mo) 

0.6897±0.0034 

1.043±0.055 

0.944±0.010 

Ri (Rq) 

0.6489±0.0013 

0.964±0.017 

0.833±0.011 

Teff, (K) 

4450±150 

5636±100 

5527±100 

Frot,l (d) 

35.1±1.0 

7.775±0.022 

20.31±0.47 

Ms (Mo) 

0.20255±0.00066 

0.362±0.013 

0.1951±0.0020 

R 2 (Rq) 

0.22623±0.00059 

0.350±0.006 

0.2150±0.0014 

Teff, (K) 


3357±100 

3226±100 

Frot,2 (d) 




[Fe/H] (dex) 

-0.20 

-0.25 

-0.34 

Closest Evolutionary Track 

Ml* (Mo) 

0.70 

1.05 

0.90 

Ms* (Mo) 

0.20 

0.35 

0.20 

Zr (Mo) 

0.010 

0.010 

0.008 

Atms.i (Gyr) 

> 13 

5.53 

10.1 

Binary Properties 

Tags (Gyr) 

2.0-3.0 

1.0-5.0 

1.5-2.5 

Fbin (d) 

41.07922±0.00007 7.44837695±0.00000021 

27.322037±0.000017 

^^bin (-^U) 

0.22431±0.00035 

0.0836±0.0014 

0.18539±0.00066 

eioin 

0.15944±0.00062 

0.0234±0.0010 

0.0510±0.0037 

Planetary Orbit 

Pp (d) 

228.776±0.029 

303.158±0.072 

240.503±0.053 

ttp (AU) 

0.7048±0,0011 

0.989±0.016 

0.7903±0.0028 

Gp 

0.0069±0.0012 

< 0.411 

0.0359±0.0088 

References 


Doyle et al. 2011; Winn et al. 2011 

Orosz et al. 2012b 

Welsh et al. 2015 


Table 1 

The KBHZ sample, Kepler binary systems harboring planets in the BHZ. Note. Quantities marked with an are model 

parameters. 


spectral distribution of each source and the response of 
planetary atmosphere to these spectra. 

To obtain the limits of the BHZ, ^Sbin given by Eq. (3), 
the equation should be averaged over a time long enough 
compared to the periods of the system. If the binary 
eccentricity is not too large, an analytical expression for 
{S) may be used: 


with n, m = 1, 2, Am = d? r‘^ and Bm = 2 d Vm- 

Finally, the BHZ limits are calculated by solving the 
equation 

Sefi,i = (5'bin)(^bin,i) (6) 

Although other authors have developed alternative 
methods to estimate the limits of the BHZ (Haghigh- 
ipour & Kaltenegger 2013; Kane & Hinkel 2013; Mason 
et al. 2013; Cuntz & Bruntz 2015; Muller & Haghighipour 
2014), all of them share the most important character¬ 
istics of the method presented here. In particular, our 
method is an improved version of the first order method 
used in Mason et al. (2013) by incorporating some key 
elements of the methods presented by Haghighipour & 
Kaltenegger (2013) and Cuntz & Bruntz (2015). An in¬ 
teresting advantage of our method with respect to others 
is the fact that for low eccentricities, the calculation of 
the BHZ limits requires finding the roots of an analytical 


formula. No numerical integrations or calculations over 
a grid of points around the system are required. This 
feature greatly improves the speed of calculation, which 
is especially important for the exploration of the vast 
parameter space of circumbinary habitability. 

We have compared the BHZ limits for all the Kepler 
binaries with known circumbinary planets calculated us¬ 
ing our method and those used in three different inves¬ 
tigations. The results are presented in Figure 1. We 
have found that the average of the relative difference be¬ 
tween the BHZ limits using the new method and that 
of Haghighipour & Kaltenegger (2013) is less than 4% 
^ The difference found by comparing our method with 
that of Mason et al. (2013) is less than 1%. The limits 
calculated by Cuntz & Bruntz (2015) are systematically 
more conservative than all of the other methods, with 
the inner limits ^10% farther out and the outer limits 
^10% closer in than the other methods. 

3.1. Evolution of Insolation and BHZ 

Generally, it is possible to assume that moderately 
separated stars in binary are formed at the same time. 
Moreover, stellar evolution shows that stars with differ¬ 
ent masses evolve at different rates and in different ways. 

® The method devised by these authors produces dynamical and 
asymmetric BHZ. Comparison with our method depends on the 
assumptions made about how to convert their asymmetrical limits 
into our symmetric ones. Part of the discrepancy may involve this 
ambiguity. 


1— — arctan 


f Bm 


(5) 
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Figure 1. Comparison of BHZ limits for all the Kepler binaries 
having circumbinary planets, calculated with four different meth¬ 
ods: Haghighipour & Kaltenegger (2013); Mason et al. (2013); 
Cuntz & Bruntz (2015) and the current paper. The limits calcu¬ 
lated with the method presented here are used as reference values 
(abcsisa values and solid black lines in upper panel). Dashed lines 
correspond to a difference of 0.1 AU (~10%). 

Thus, in order to investigate habitability around a bi¬ 
nary, stellar evolution models for each star must be used 
to determine the time dependence of the BHZ. 

It is important to ask whether the fundamental proper¬ 
ties of low-mass stars, such as effective temperature and 
luminosity used to calculate the BHZ, evolve in binaries 
exactly as they do in single-stars. For moderately sep¬ 
arated binaries, (periods larger than 8 days, stellar sep¬ 
arations larger than tens of stellar radii), such as those 
studied in our sample, the Roche lobes of the stars are 
one order of magnitude larger than the stellar diameters. 
For instance, in the case of the tightest binary in our 
sample (Kepler-47), the filling factor of the primary (ra¬ 
tio of stellar radius to Roche lobe size) is ^0.2 during 
PMS (when the star is largest) and ^0.1 during main- 
sequence (Eggleton 1983). At those filling factors the 
shape of the star is only slightly modified by the com¬ 
panion and its interior and atmospheric properties are 
practically the same as for a single-star. 

HZ evolution is also a natural outcome of the evolu¬ 
tion of single-stars. In the same way that the continuous 
HZ is defined for single-stars (Kasting et al. 1993), it 
is generalized to binaries and a continuous circumbinary 
habitable zone (CBHZ) is defined. 

In Figure 2, the evolution of the BHZ for Kepler-47 and 
its associated CBHZ is shown. We use stellar evolution¬ 
ary models, at each age, to calculate the instantaneous 
properties of both stars. Then the limits of the HZ at 
that time are determined. Shown here are the the Recent 
Venus and Early Mars criteria. The outer CBHZ limit is 
defined by the instantaneous outer edge of the BHZ as 
calculated at the time when the binary flux is minimum. 



Figure 2. Kepler-47 is shown as an example of the evolution 
of the BHZ (green area). BHZ limits are shown as a function of 
age (solid blue and red lines). The HZ for a single-star with a 
mass equal to the primary is also shown (dashed lines). Since the 
system is composed of a solar-like star and an M dwarf the limits 
for the single primary and the binary are almost the same. The 
CBHZ is highlighted (shaded gray region; see text for explanation). 
The semimajor axis of Kepler-47c, is shown as a horizontal line. 
Neglecting eccentricity effects, the habitability lifetime of this orbit 
is roughly 3.5 Gyr. 


We have tested that this time is close to the time of zero- 
age main-sequence (ZAMS) of the primary star, which by 
definition is the more massive component of the system. 
Planets in the CBHZ remain inside the BHZ during the 
entire main-sequence stage of the primary as it is the star 
with the shortest lifetime. The definition of the CBHZ 
inner edge is trickier. We define it as the position of the 
inner edge of the instantaneous BHZ calculated when the 
primary star abandons the main-sequence (close to the 
end of the nuclear hydrogen-burning phase). 

In order to assess the present insolation conditions in 
the KBHZ sample we plot, in Eigure 3, the instantaneous 
BHZ calculated in the middle of the range of estimated 
age (See Table 1 for ages and references). Along with 
BHZ limits, the insolation and photosynthetic photon 
fiux density (PPED) are shown as a function of planetary 
orbital phase. See Mason et al. (2015) for details on 
PPED calculations. 

In the BHZ diagrams of Eigure 3, the planetary or¬ 
bit has been placed among different landmark distances: 
BHZ and CBHZ edges, binary orbit, equivalent Earth 
insolation distance, and last but not least, the critical 
stability distance Ucrit, defined as the minimum distance 
where test particles have stable circumbinary orbits. We 
compute this critical distance with the semi-empirical re¬ 
lationship of Holman & Wiegert (1999): 

(^crit ~ (1-60 + 5.10ebin “ 2.22ebin + (7) 

4.12ya — 4.276bin/^ — 5.09/i^ T 
4.61ebin/^ ) ^bin* 
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^in,Rv=0-36 AU, /i„,RG=0.46 AU, U,mg=0.83 AU, /out,EM=0-88 AU 
a,rit=0.65 AU, /i„,,o„t=0-65 AU, U,cont=0.85 AU 



0-50 AU Kepler-16 




orbital phase 




Figure 3. BHZs, insolation and PPFD (see Mason et al. (2015) for an explanation) for the KBHZ sample. Left: Along with the BHZ 
instantaneous limits, the binary orbit, the critical distance (dashed black line), the planet orbit (solid black line), the limits of the continuous 
BHZ (orange solid lines), and distances of binary orbit resonances (concentric gray circles) are shown. The distance at which the equivalent 
Earth insolation is received is marked with a dotted black line. Right: Insolation and PPFD are calculated in units of present Earth levels 
not only at the distance of the planet (blue and red solid lines, respectively) but also at all locations within the BHZ (shaded areas). 
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Here fi = M 2 /{Mi + M 2 ), and abin and ebin are the 
binary semimajor axis and eccentricity, respectively. Al¬ 
though dynamical instability likely occurs outside the 
critical limit as well, as also discussed by Holman & 
Wiegert (1999), all the planets in our sample are safely 
beyond this critical distance. However, in Kepler-16, Ucrit 
lies inside the BHZ, restricting the instantaneous exten¬ 
sion of the BHZ as well as the CBHZ; see top left panel 
of of Figure 3. When Ucrit is larger than the inner edge 
of the continuous BHZ we use its value as the inner edge 
distance. 

In the instantaneous BHZ diagrams of Figure 3, dis¬ 
tances at which orbital resonances with the binary oc¬ 
cur are also shown. It is very interesting to notice, and 
likely not coincidental, that in both cases, Kepler-16b 
and Kepler-453b, the planet remains entirely in between 
two resonances throughout its orbit. This supports the 
idea that orbital binary resonances may play a role in 
the formation and final architecture of circumbinary sys¬ 
tems. Notice also that Kepler-47c crosses many tightly 
spaced resonances, not plotted in its BHZ diagram in 
Figure 3 for clarity. 

The planet Kepler-16b has a nearly circular orbit just 
inside the maximum greenhouse limit. Kepler-47c has 
a highly eccentric orbit and leaves the BHZ for part of 
its orbit. It does, however, currently spend most of its 
time within the HZ, especially given the slower velocity 
at apastron (see insolation diagram in the middle right 
panel of Figure 3). The planet Kepler-453b has a slightly 
eccentric orbit and always remains in the HZ, just inside 
the runaway greenhouse limit. 

According to our calculations, of all the circumbinary 
Kepler planets, only Kepler-16b is in the CBHZ, i.e. it 
will he inside the BHZ during the whole lifetime of the 
primary star, which is in turn larger than the age of the 
universe. This condition would seem at first very fa¬ 
vorable for the emergence and evolution of complex life. 
However other factors could limit the long-term habit¬ 
ability within environments surrounding Kepler-16 and 
similar binaries (see Section 8). 

Kepler-47c and Kepler-453b do not remain inside the 
BHZ for the main-sequence lifetime of the primaries. 
However, both remain inside (although very close to 
the inner edge of) the BHZ for astrobiologically relevant 
times (3.5 Gyr and 4 Gyr respectively). If we believe 
stellar evolution models, Kepler-47c is about to leave or 
has just left the most optimistic BHZ, being currently in 
a Venus-like state in terms of stellar insolation. 

In order to calculate BHZ evolution and other key 
quantities for models, the stellar properties available in 
the PARSEC vl.2 evolutionary tracks are used (Bressan 
et al. 2012, 2013) In addition to providing good cov¬ 
erage in metallicity and mass, the PARSEC tracks sample 
the PMS stages of stellar evolution well (and of course 
also the main sequence), a key feature when attempting 
to study the early evolution of stellar rotation. More¬ 
over, these models reflect improvements for low-mass 
stars (Chen et al. 2014). In the BHMcalc tool we pro¬ 
vide access to other evolutionary model results. 

4. THE EVOLUTION OF ATMOSPHERES AS A KEY 
FACTOR FOR ASSESSING HABITABILITY 

http://people.sissa.it/ bressan/parsec.html 


So far we have calculated insolation as a way to assess 
habitability conditions in circumbinary planets. But in¬ 
solation is not the only factor constraining the capability 
of a planet to sustain life (for a review of many other 
factors see Kasting 2009 and references therein). Among 
other endogenous and probably exogenous factors, the 
presence of a properly dense atmosphere with the right 
amount of water is mandatory. 

The solar system is a good example of a planetary sys¬ 
tem where among three planets (and a sizeable moon), 
located close or inside the HZ, two of them, Venus and 
Mars (plus the Moon), are un-inhabitable mainly due 
to factors related to the composition and mass of their 
atmospheres. 

The evolution of a planetary atmosphere is driven by 
the influx of high-energy radiation and plasma from the 
host star. This has been a matter of research in the last 
two decades (Kasting 1996; hammer et al. 2007, p. 127; 
hammer et al. 2009, p. 399; hammer et al. 2010; ham¬ 
mer 2013; Tian et al. 2013, p.567). Although it is not yet 
clear which factors determine the survival of Earth-like 
atmospheres or drive their composition to favor habitable 
conditions, multiple lines of evidence show that the in¬ 
teraction of a planet with its young host star is a strong 
factor in determining the fate of potentially habitable 
planet atmospheres (hammer et al. 2007, p. 127; ham¬ 
mer et al. 2009, p. 399; Tian 2009; Zendejas et al. 2010; 
Zuluaga et al. 2013). 

Two key factors are involved in this interaction. The 
first is the level of X-rays and extreme-ultraviolet (EUV) 
radiation, collectively called XUV radiation, incident on 
the planet. These photons heat up the upper atmosphere 
and under certain conditions could lead to massive loss of 
volatiles (Tian 2009; Tian et al. 2013, p. 567). The flux 
of charged particles, mainly protons, and their associated 
magnetic fields (SW) is the other factor. The direct in¬ 
teraction between this continuous flow of particles and 
fields and planetary atmospheres strongly impacts un- 
magnetized planets and moons (e.g., Venus, Mars and 
the moon). Resulting nonthermal processes may remove 
large fractions of planetary atmospheres, or in extreme 
cases they can strip them almost completely (Zendejas 
et al. 2010). Mars is the best known example of the ef¬ 
fect of this kind of aggression eliminating habitability. 
Recent data from the MAFEA Mission (Rahmati et al. 
2014; Jakosky et al. 2015) are showing the increased level 
of importance that SW interaction with the atmosphere 
has in shaping the evolution of Mars’s capability to have 
liquid water on its surface (habitability). 

Modeling the interaction between the radiation and 
plasma fields of the star and, in this case, the binary 
system is key in assessing the habitability of all plan¬ 
ets. However, since binaries are relatively novel environ¬ 
ments, at least for planetary systems, we need to under¬ 
stand and learn to model aspects of stellar evolution that 
have been applied for years to single-stars. In the follow¬ 
ing sections we develop a theoretical framework to model 
the evolution of stellar activity and its role in the chang¬ 
ing conditions faced by the atmosphere of a circumbinary 
habitable planet. 

5. STELLAR ROTATION AND MAGNETIC ACTIVITY IN 
BINARIES 

5.1. Evolution of Stellar Rotation 


n 


The evolution of stellar rotation has been a matter of 
intense research over the last several decades, not only 
in single but also in binary stars (for a recent review see 
Bouvier 2013 and references there in; for the evolution 
of rotation in binaries see Claret et al. 2005). 

The large diversity of stellar rotation periods observed 
in stars of all ages (especially PMS stars) has driven the 
development of a diverse array of semi-empirical and phe¬ 
nomenological models describing the evolution of stellar 
AM (Kawaler 1988; Chaboyer et al. 1995; Matt et al. 
2012). These models have been extensively tested against 
samples of stars (both single and binary) in young stel¬ 
lar clusters and more recently in low-mass field stars (Ir¬ 
win et al. 2011). The details of these models depend on 
weakly constrained parameters; however, a general con¬ 
sensus has been reached about the general forces and 
timescales driving stellar rotation evolution. 

Here we apply and extend some of these models for 
stars in moderately separated binaries (i.e. binaries with 
similar masses and orbits to those of the KBHZ bina¬ 
ries). For details concerning the physical motivations, 
history of development, and success of these models at 
reproducing observations, please refer to Bouvier (2013) 
and references therein. 

Although the applicability of these models to stars in 
binary systems has not been systematically tested, we 
assume that some of the underlying mechanisms will op¬ 
erate in the stellar components of these binaries as if 
they were isolated stars. A detailed test of this hypothe¬ 
sis should be pursued in future research efforts. Still, for 
the range of orbital separations and stellar masses con¬ 
sidered here, the evolution in insolation hypothesis re¬ 
main as a good one for the same arguments given before 
to support the applicability of single-star evolutionary 
models. 

In general, the AM of the ith layer of a star evolves 
while obeying Newton’s angular second law: 


Here Ji = is the instantaneous AM of the layer, 
li its moment of inertia, its instantaneous angular 
velocity, and 7i is the net torque over the layer. 

In the case of low-mass stars (AA < I. 3 M 0 ), it is cus¬ 
tomary to assume that the star has two distinct interior 
layers: a radiative inner core {i = core) and a convective 
envelope (i = conv). It is also assumed that both layers 
rotate independently as rigid bodies (i.e. no differential 
rotation). 

The convective envelope is subject to two different ex¬ 
ogenic effects: (1) tidal breaking, Ttid, and (2) loss of AM 
via a magnetized SW Tsw- 

To calculate Ttid, we use the formalism developed by 
Hut (1981). Accordingly, a target star of mass Mtarg and 
radius i^targ, rotating at an instantaneous rate Q. and 
affected by a secondary star of mass Mfieid ia an orbit 
with semimajor axis a, eccentricity e, and mean angular 
velocity n, feels an effective tidal torque given by 


Ttid 


3^2-^ f ATfield \ ^ f Rtarg \ ^ 

^^^diss \ Aftarg J \ ^ J 


(9) 


^ (1 - e2)« 


hie’) - (1 - e’'fl’hie’)^ 


Here k 2 is the apsidal motion constant that quantifies 
the response of the star to tidal distortion, and it de¬ 
pends in general on the internal mass distribution of the 
star, 

fdiss (Mtargi^targ/^targ)^/^ IS the timCSCalc 
for dissipation of tidal energy inside the target star, and 

K, = y^//(A/targ^targ) ^^Lc gyration radius. Eccentric¬ 
ity functions / 2 (e^) and / 5 (e^) are defined by Eqs. (2) 
and (3) in Hut (1981). 

We assume that while tides are raised over the whole 
star, only tides within the convective layer contribute to 
AM dissipation (Zahn 2008). As a result, the torque in 
Eq. 10 only affects the envelope of the star. Alterna¬ 
tively, we can assume that tides raised in the radiative 
core are dissipated on timescales much larger than that 
in the convective zone, rendering the resulting torque 
comparatively small. 

The parameter /diss quantifies the efficiency of tidal 
energy dissipation. Although a fundamental theory of 
this process capable of correctly explaining observations 
of close binaries is still lacking (see Claret 2011 for a 
discussion) we assume for simplicity that /diss ~ 1, as is 
used, for example, in Claret (2012). Alternative theories, 
predict values of /diss within one order of magnitude. 
Thus, for instance, the turbulent dissipation theory of 
Zahn (2008) predicts a value /diss ~ 3.5. 

Eor nonnegligible eccentricities, tidal braking could 
drive stars to a pseudosynchronous final state where ro¬ 
tational angular velocity is related to, but not exactly, 
the average orbital angular velocity (Hut 1981), 


_ ^sync _ 1 + 

- (I_e2)i(l + 3e2 + |e4) 

Eor binary eccentricities e in the range 0 — 0.5, 1 < 
Usync <2.8, implying that provided enough time, stellar 
components of a binary will be locked in a rotational 
state where their periods are close to that of the binary 
period. This result will not depend on the mass loss 
history of each star. Once locked, the magnetic activity 
of the stars becomes frozen at a given value, which in 
some cases could mimic the behavior of a star that is 
much younger or older (Mason et al. 2013). 

The second exogenic effect removing AM from the con¬ 
vective layer is mass loss via magnetized SWs. Not only 
does the mass leaving the surface of the star carry AM, 
but in magnetically active stars the flow of charged par¬ 
ticles is also coupled with the corotating magnetic field, 
up to few to several tens of stellar radii. This leads to a 
loss of AM that is several orders of magnitude larger than 
that in stars without a convective envelope (Schatzman 
1962). 

Although a successful, comprehensive, and self- 
consistent model of magnetized SW AM loss has not yet 
been developed, a semi-empirical formula, developed by 
Kawaler (1988) and modified by Chaboyer et al. (1995), 
succeeds at producing results consistent with observa¬ 
tions for a large range of stellar masses and ages. Accord¬ 
ingly, the torque exerted by the magnetized SW scales 
with angular velocity 11, stellar mass M^, and radius 


X 
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following 


/ D \ 1/2 / IWf \ “1/^ 

Tsw = ^ min(f7, f^sat)^ (/t^y 

( 11 ) 

Kc and llgat are two free parameters that can be 
adjusted to reproduce the scale of observed rotational 
rates and their distribution, respectively. In order to re¬ 
duce the number of free parameters in our comprehensive 
model, flsat is calculated following the scaling relation¬ 
ship 


= ( 12 ) 

\ '^conv J 

where Tconv is the convective overturn time (see Section 
5.2.2). This scaling relationship has proven to be useful 
in reproducing the distribution of rotation periods for 
low-mass field stars (Irwin et al. 2011). 

More recently, success in modeling the rotational rate 
distribution of stars in young clusters was achieved (Matt 
et al. 2012; Gallet & Bouvier 2013) using a semianalytical 
fit of MHD simulations to the original model of Kawaler 
(1988) along with a recent model for stellar mass loss 
(Cranmer & Saar 2011). However, its application is still 
restricted to stars with mass nearly identical to the Sun. 
Our interest here is in stars ranging from 0.2 to 1.05 Mq; 
see Table 1. A comparison between the results obtained 
with the semi-empirical formula used here (Eq. 11) and 
those obtained with the formula by Matt et al. (2012) is 
left for a future work, but it is not expected to alter the 
major conclusions of the current work. 

As the convective envelope loses AM, a difference in 
angular velocity between the envelope and the radia¬ 
tive core is established. A variety of physical processes, 
such as magnetic fields, hydro dynamical instabilities, 
and gravity waves, may redistribute AM in the stellar in¬ 
terior, reducing this difference with time. To model AM 
redistribution we follow, as is customary, the prescrip¬ 
tion by MacGregor & Brenner (1991), by estimating the 
mutual torque due to the difference in rotation rates as 
given by 


|7aj| =- (13) 

'^ce 

The AM exchanged is AJ = (/conv-^core — 
IcoreJconv)/hot^ and Tee is a free parameter called the 
core-envelope coupling timescale. 

In summary, the AM of the convective envelope and the 
radiative core of each star in the binary evolve following 
the equations 


dJ 


dt 

dJrad 


dt 


Ttid + Tsw + Taj 
—TaJ‘ 


(14) 

(15) 


The moment of inertia of both the convective envelope 
and the radiative core evolves in time and affects the 
rotational history of the whole star. PMS evolution is 
especially important in this respect. During PMS, stars 
contract the most. A pragmatic approach to account for 


this early contraction is to integrate the angular velocity 
instead of the AM using 


dQi 1 dJi ^ dli 
dt li dt * dt 


(16) 


where the term —Qi{dli/dt) accounts for the contraction 
(expansion) of the respective layer. It is important to 
stress that although the larger changes in li happen dur¬ 
ing the PMS, we are also integrating Eq. (16) during the 
main sequence. 

Erom a purely dynamical point of view this contrac¬ 
tion term is associated with two discrepancies. The rapid 
contraction of the star during the first stages of stellar 
evolution is compounded with the gain of specific AM by 
accreting gas. Angular acceleration of the star may con¬ 
tinue even to near breakup velocities (Prot ^0.1 days). 
Observations, however, show that the rotation rates of 
PMS stars are orders of magnitude slower (with periods 
of 2 — 10 days). Although still debated, the solution to 
this discrepancy may rest on the assumption that stars 
transfer this AM excess to their accretion disk. This disk 
locking persists for a timescale Tdisk ^2 — 5 Myr, ending 
with disk dissipation. 

The second discrepancy arises when considering the 
value and sign of the contraction term for the radiative 
core. In stellar evolution simulations, the core suddenly 
forms at around 10 — 100 Myr and grows to its final size 
over a time an order of magnitude smaller. As a result, 
the moment of inertia of the core increases rapidly, po¬ 
tentially producing a sudden rotational breaking of this 
layer. Even after considering the exchange of AM during 
the growth of the core, as is done in Allain (1998), this 
sudden braking still occurs. Although unobservable, a 
sudden deceleration of the radiative core would produce 
a signature in the evolution of the envelope and would 
create a discrepancy with observations. 

To be consistent with previous results, we assume that 
the contraction term for the core in Eq. 16 is the same as 
that of the convective envelope. An investigation of how 
this apparent deceleration in the radiative core seems 
to be prevented is a critical matter in need of further 
research. 

Eor each system studied here, numerical integrations 
of Eqs. (16) were performed for both stars in the bi¬ 
nary including their mutual tidal interaction. Eor that 
purpose we take as inputs the values of key stellar struc¬ 
ture properties as a function of time (e.g. apsidal motion 
constants, moments of inertia, and convective zone thick¬ 
ness), from stellar evolutionary models (Bressan et al. 
2012, 2013) as well as from interior structure models 
(Baraffe, 2014, personal communication). 


5.2. Magnetic Activity 

The ultimate purpose of calculating instantaneous stel¬ 
lar rotational periods is to calculate the magnetic activ¬ 
ity as a function of time. Magnetic phenomena in the 
atmospheres of cool stars are driven by the action of a 
dynamo in their convective zones that in turn is respon¬ 
sible for producing SWs and high-energy radiation (XUV 
emission). 

Progress in understanding and characterizing these 
processes in the Sun, allowed S. Granmer and S. Saar 
to develop a self-consistent model of SW acceleration 
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applicable to other cool stars (Cranmer & Saar 2011). 
Their model accurately predicts stellar activity observ¬ 
ables such as the average photospheric magnetic field 
strength and mass loss rates as a function of several basic 
stellar properties such as mass, radius, luminosity, and 
rotational period. 

Their publicly available routine BOREAS is applied here 
in order to calculate the instantaneous mass loss rate of 
each star in the binary system. For details on the input 
physics, assumptions, and free parameters used by the 
models of Cranmer and Saar, refer to Cranmer & Saar 
( 2011 ). 


5.2.1. SW 

Armed with the mass loss rate M, the properties of 
the SW are calculated in a straightforward manner. We 
assume that the wind reaches its terminal velocity usw 
at a distance from the stellar center, being equal to 
the local escape velocity: 


'^sw 


2GM^ 


n 


(17) 


A distance rt = 2.5i?* is assumed in order to reproduce 
the current average solar wind condition as observed at 
Earth distance. 

By mass conservation, the particle mean density as 
measured at a distance r from the star is given by 


M . , 

where fi is the mean molecular weight of the wind parti¬ 
cles and Trip is the proton mass. For a fully magnetized 
SW li = 0.6 (Cranmer 2008; Cranmer & Saar 2011). 

This stream of charged particles continuously impacts 
the atmosphere of planets, at a rate given by 


Fsw = ^sw'^sw- (19) 

At the distance of the HZ, the wind is supersonic and 
has a negligible thermal pressure. As a result, it exerts 
a dynamic pressure over planetary magnetospheres given 
by 


ftw = M^p^sw('^sw + '^orb) (20) 

where Uorb is the orbital velocity of the planet. 

The Cranmer & Saar model has been developed in or¬ 
der to describe activity in single-stars. The previous 
assumptions concerning how to convert mass-loss rates 
into SW properties, may work well in describing SWs 
for moderately spaced binaries, as well as isolated stars, 
but the applicability of these models to binaries has not 
been tested. Additionally, other effects could arise, mod¬ 
ifying the underlying assumptions on which this model 
is built. For example, a star could illuminate the atmo¬ 
sphere of the companion, thereby modifying the vertical 
atmospheric structure. Coronal X-ray heating may play 
a critical role in some cases. The acceleration due to 
gravity is also modified by the presence of a close com¬ 
panion. Stellar wind shocks could also arise, modifying 
the simple model described herein. We assume that the 
acceleration mechanism working in single stars provides 


a valid constraint for stars in the moderately separated 
binaries under consideration. 

At distances large enough from the binary or in very 
disparate binaries {q ^ 1), it is reasonable to assume 
that the SW flux and dynamic pressure are simply the 
sum of the quantities produced by each star, i.e. 


^SW,bin " 

- ^SW,1 + ^SW,2 

(21) 

^SW,bin - 

" ^SW,1 + ^SW,2- 

(22) 


For close orbiting planets, these simple assumptions 
are expected to break down under realistic conditions in 
some binaries. Recently, the conditions under which the 
interaction between the winds of binaries could strongly 
interact, rendering this approximation imprecise, have 
been studied using MHD simulations. The results of 
these simulations suggest that even in an extreme case 
of solar-mass twins, the SW flow is only strongly mod¬ 
ified at close circumbinary distances. Within the BHZ 
the simple approximation given in Eq. (21) is applicable 
for the purposes of studying the plasma environment of 
potentially habitable planets. 

5.2.2. XUV Emission 

One of the most detrimental effects that planets face 
is early stellar activity. During the early history of plan¬ 
ets, the high-energy radiation flux far exceeds that which 
they receive later on. In the case of Earth for example, it 
has been estimated (Guinan & Engle 2009) that 4.5 Gyr 
ago the X-rays flux was 10-100 times larger than present 
Earth level (hereafter PEL) owing to solar coronal activ¬ 
ity. 

The effect that such flux has over time on planetary 
atmospheres ranges from simply a modification of the 
upper atmospheric chemistry (Tian et al. 2008; Segura 
et al. 2010; Hu et al. 2012, 2013; Hu & Seager 2014) to 
rapid photoevaporation via induced heating of the exo¬ 
sphere (Lammer et al. 2003; Lammer et al. 2009, p. 399; 
Lammer et al. 2010; Zendejas et al. 2010). Both phe¬ 
nomena could have very detrimental effects on planetary 
habitability. Estimating the flux of high-energy radiation 
is thus an essential element in the assessment of habit¬ 
ability of any planetary environment. 

The emission of X-Ray and EUV radiation (1-91.2 nm) 
and their dependence on stellar age and rotational period 
have been extensively studied in the case of single-stars. 
A strong anticorrelation between age and XUV luminos¬ 
ity has been observed (Ribas et al. 2005; Penz & Micela 
2008; Penz et al. 2008). 

Binary stars follow a different set of rules. First, ro¬ 
tation and age are not necessarily correlated, especially 
when the binary period is short (of the order of days) 
or if the binary eccentricity is larger than about 0.1, ow¬ 
ing to tidal torques. Second, and as explained before, 
strong magnetic interaction between the stellar compo¬ 
nents in the first tens of Myr could drive rotational evolu¬ 
tion and hence activity toward values that are uncommon 
in single-stars. 

In spite of the most common recipes used to estimate 
XUV levels as a function of stellar age (see, e.g., Zuluaga 
et al. 2012), which will not be applicable in the case 
of binaries for the reasons given above, we use here an 
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approach consistent with the methods used to calculate 
mass loss rates, i.e. using the activity proxies provided 
by the models of Cranmer & Saar (2011). 

Using a catalog of ^1600 stars whose rotation and X- 
ray brightness have been measured, Vardavas (2005) and 
more recently Wright et al. (2013) showed that X-ray 
luminosity is strongly correlated with the Rossby number 
(the ratio of inertial to the Coriolis force in the convection 
zone). The Rossby number is calculated as 

Ro = Prot/Tconv (23) 

Here P^ot is the rotational period (which determines 
Coriolis forces) and Tconv (the convective overturn time) 
is the typical timescale for convective motions (inertial 
forces). 

Wright et al. (2013) found the following empirical re¬ 
lationship: 


Rx 


Lx 

Lho\ 


= Px,sat X 


1.0 if Ro < ROsat 

(Ro/Rosat)~^‘^ otherwise 


(24) 

Here Rogat and Px,sat are fitting constants quantifying 
the Rossby number at which the emission is saturated 
and the constant ratio of X-ray to bolometric luminos¬ 
ity observed at saturation levels. Using their catalog, 
Wright et al. (2013) found that (Rosat,Px,sat) are be¬ 
tween a high-level range of (0.17,0.001075) and a low- 
level range of (0.09, 0.00051). 

It is interesting to notice that at a given value of Ro the 
X-ray luminosity could be predicted using the empirical 
law in Eq. (24) with an uncertainty of a factor of 2 
even among a large range of ages, stellar masses, and 
rotational periods. 

In order to incorporate this result into our model, we 
need to predict the Rossby number of a star at a given 
stage of its rotational evolution. For that purpose, the 
routines of Cranmer & Saar (2011) that provide the con¬ 
vective overturn times as a function of stellar effective 
temperature are used. 


6. THE SOLAR REFERENCE MODEL 


In order to calibrate binary star models, we have calcu¬ 
lated the evolution of rotation and activity for an isolated 
solar-mass star. The results are presented in Figure 4. 
We call this the Solar Reference Model (SRM). 

The free parameters of the model that were described 
in previous sections have been chosen to reproduce the 
observed properties of the Sun, namely, the period of 
rotation, mass loss rate, and XUV luminosity at present 
times. 

Since activity and rotation of the Sun during the early 
phases of solar system evolution are poorly constrained, 
three scenarios are considered: (1) a nominal scenario 
(initial rotation period Prot,ini = 7 days. Tee = 7 Myr), 
(2) an initially fast-rotating Sun (Prot,ini = 1-5 days. 
Tee = 15 Myr) and (3) an initially slowly rotating sun 
(Prot,ini = 12 days. Tee = 1 Myr). In each case we use a 
disk locking time Tdisk = 2 Myr, in agreement with the 
estimations of the solar nebula dissipation time (Krot 
et al. 2005). For each scenario the value of the braking 
constant Kc has been set so as to reproduce the prop¬ 
erties of the Sun at r = 4.56 Gyrs (see legend in upper 
panel of Figure 4). 



Figure 4. Rotation and activity evolution of a solar-mass star in 
three different scenarios: a slow, nominal, and fast initially rotat¬ 
ing Sun (see text). Upper panel: rotational velocity as a function 
of time. The measured rotations of five solar-like stars (masses 
between 1.01 and 1.10) and the Sun itself are included for refer¬ 
ence (the periods of rotation of the reference stars are taken from 
Ribas et al. 2005). Middle panel: XUV luminosity as a function 
of time in terms of the present-day solar X-ray luminosity. Lower 
panel: evolution of mass-loss rate. 


Using rotation periods, the mass-loss rate, calculated 
from the BOREAS routine, as a function of time for each 
scenario is shown in the lower panel of Figure 4. The 
XUV luminosity has been calculated in each case us¬ 
ing Eq. 24 and by assuming a critical Rossby num¬ 
ber Rosat = 0.16 and a saturation X-ray luminosity of 
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Figure 5. Integrated XUV (upper panel) and SW (lower panel) 
fluxes as measured at Venus, Earth and Mars orbits. Integration 
is performed starting at 10 Myrs for the primary atmospheres and 
700 Myrs (inset panels) for secondary atmospheres. 

^x,sat = 7.4 X 10“^. These values are chosen in order to 
reproduce the observed X-ray luminosity of the Sun. 

With the calculated XUV luminosities and mass loss 
rates in hand, the SW and XUV fluxes measured at 
the distance of Venus, Earth, and Mars are determined. 
From these fluxes, we calculate the integrated effect on 
a planetary atmosphere: 

^xuv,sw('r; Tini) = J Txuv,sw(^)<^^ (25) 

Integrated values of XUV and SW fluxes are shown 
in Figure 5. Integrated fluxes rapidly reach a maximum 
around 200 — 400 Myrs in the case of SW and ^1 Gyr 
in the case of XUV radiation. At about those times, 
any hypothetical primitive planetary atmosphere (pro¬ 
toatmosphere) has received most of the incoming flux 
that it will ever receive. 

If we assume that a “secondary” atmosphere is de¬ 
gassed or left over from the formation process (see, e.g., 
hammer 2013 and hammer et al. 2013 and references 
therein), we can also calculate the cumulative effect 
that SW and XUV radiation has on these secondary en¬ 
velopes. In the inset panels of Figure 5, we plot the inte¬ 
grated flux starting at 700 Myrs (around 3.8 Gyrs before 
present) and up to 8 Gyrs. Again, most of the integrated 


effects reach an asymptotic value after 200 — 300 Myrs 
(around 1 Gyr after formation). 

Values of the integrated particle and photon fluxes can 
be converted into total atmospheric mass loss from un¬ 
magnetized planets. However, converting one quantity 
into the other is nontrivial given the complexity of the 
nonthermal and thermal mass loss processes involved. 
For present purposes, first-order estimations can be ob¬ 
tained with the phenomenological formalisms of Zende- 
jas et al. (2010) (for atmospheric stripping) and Sanz- 
Forcada et al. (2011) (for XUV-driven erosion). 

Using the asymptotic values of the integrated fluxes 
in Figure 5, a total solar-driven mass loss of the order of 
0.1 —0.2 bars for a hypothetical “primary” atmosphere of 
Mars is obtained. Venus could have lost 1.5 — 3 bars from 
any protoatmosphere as a result of the same mechanism. 
In both cases, and to be conservative, we assumed that 
both planets started with G 02 -rich atmospheres^ 

For secondary atmospheres (i.e. any atmosphere re¬ 
maining or degassed after several hundreds of Myr), we 
estimate that Mars and Venus were subject to mass losses 
of the order of 0.06 — 0.1 bars and 0.8 — 1.5 bars, respec¬ 
tively. These numbers are in agreement with previous 
estimations (see, e.g., Kulikov et al. 2006). 

In the case of Mars, which currently has a 6 mbar at¬ 
mosphere, an estimated mass loss of ^100 mbar is fairly 
compatible with an early lost of “habitable” conditions 
(the presence of surface water). Venus, having a present 
^100 bar atmosphere, seems to have barely lost 1% of 
its mass via direct stripping from solar wind. 

If, on the other hand, we assume that terrestrial plan¬ 
ets started with massive hydrogen/helium envelopes (see 
hammer 2013 and references therein), subject to photo¬ 
evaporation from the absorption of XUV radiation, we 
estimate (from the integrated XUV fluxes) that the pro- 
toatmosphere of Venus could have lost 0.3-1% of its total 
mass in the first 300 Myrs. This is of the same order as 
the mass a planet with its size could have accumulated 
from the planetary nebula or from other early outgassing 
processes (Elkins-Tanton & Seager 2008). 

In summary, and recognizing that other complex en¬ 
dogenous and exogenous factors are involved in the evo¬ 
lution of the atmosphere of terrestrial planets and moons 
(vulcanism, degassing, impacts, etc.), integrated SW and 
XUV fluxes calculated for the reference solar model as 
well as for the KBHZ sample^ are fairly compatible with 
the state of affairs in the solar system. In consequence, 
our model can be used to place first-order constraints on 
the critical level of stellar aggression that planets around 
single and binary stars could potentially withstand, be¬ 
fore losing their water inventory and becoming uninhab¬ 
itable. 

7. ROTATIONAL EVOLUTION OF KEPLER 
CIRCUMBINARY PLANETS 

The model developed and tested in previous sections is 
applied in order to calculate the evolution of rotation and 

® The terms “primary” atmosphere, “protoatmospheres,” and 
“secondary” atmospheres should not mislead the main aim of this 
work, which is to estimate the flux of particles and XUV radiation 
experiencied by planets in the early phases of stellar evolution. For 
a thorough discussion regarding the mass, composition, and evolu¬ 
tionary phases of primitive atmospheres around Earth-like planets, 
please refer to Lammer (2013) and references therein. 
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activity of the KB HZ sample. Figure 6 shows rotation 
rate as a function of time for each star in the three binary 
systems. 

In each case, the equation of motion (Eq. 16) including 
the effect of tidal interaction is integrated. For primary 
stars, we use the same parameters of the rotational evolu¬ 
tion model (flsat 5 Pini) as those of the solar reference 
nominal model. Using the same parameters for the sec¬ 
ondaries (all of them are low-mass stars M < O.SGMq) 
produces anomalously low rotational rates at late times. 
In order to be consistent with observed periods or ro¬ 
tation of field M dwarfs, we use a much lower braking 
parameter, Kc = 10^^ for companions. With this ad¬ 
justment the rotation periods of Proxima Centauri and 
Barnard’s Star are both accurately predicted. 

Although there are discrepancies between the rotation 
periods of the primary star measured photometrically 
and spectroscopically (blue and cyan error bars in Fig¬ 
ure 6, respectively), the model (blue line) predicts the 
observed values reasonably well considering the large un¬ 
certainty box. 

The observed period of Kepler-16A, P^ot = 35.1(35.8) 
days, is very close to that expected if the star is syn¬ 
chronized with the binary orbit. For the present ec¬ 
centricity of the system the pseudosynchronous period 
is Psync = 35.6 days (see Eq. 10). However, the tidal 
torque, even for the shortest expected value of tidal dis¬ 
sipation, is too low to achieve synchronization during the 
first couple of Gyrs (Nntice that the rotational evolu¬ 
tion calculated without tidal interaction is shown, using 
dashed lines, in Figure 6), for reference. Interestingly, 
the rotation evolution model also predicts a rotation pe¬ 
riod very close to the pseudosynchronous period for the 
estimated median age. This may be considered a major 
accomplishment for the proposed model. 

The case of Kepler-47A is also very interesting. The 
period of rotation measured photometrically, Prot,photo = 
7.7 days, is close although larger than the pseudo- 
syncronization period of Psync = 7.4 days. Indepen¬ 
dently, a spectroscopically determined velocity of rota¬ 
tion predicts a slightly larger rotation period Prot,vei = 
11.9 days (provided that the inclination angle is close to 
90°). Stars in this system are close enough to have their 
rotation periods substantially affected by tides. Our ro¬ 
tation evolution model predicts that within the estimated 
range of ages for the system (1 — 5 Gyr), the rotation pe¬ 
riod will reach a maximum of ^10 days. This value is 
in between the independently measured periods. The 
primary star in the model is not synchronized as early 
as first-order models suggest. Instead, once all the ef¬ 
fects modifying the rotational rate during the PMS are 
taken into account, tides become relevant only after sev¬ 
eral hundreds of Myr. 

The case of Kepler-453 is more conventional, with tides 
having a negligible effect on rotation during the relevant 
stellar evolution stages. The nominal model predicts a 
period of rotation slightly larger than the photometri¬ 
cally measured one, but nearly identical to the spectro¬ 
scopically measured period. These three cases provide 
evidence that rotational evolution is complex yet pre¬ 
dictable by the current model. 

8. RADIATION AND PLASMA ENVIRONMENT OF 
KEPLER CIRCUMBINARY PLANETS 


Kepler-16 



10 '^ 10 '^ 10 '^ 10 ° 10 ^ 


r (Gyr) 
Kepler-47 



Kepler-453 



Figure 6. Rotational evolution of each star in the Kepler bi¬ 
naries with an HZ planet, KB HZ sample. Blue lines correspond 
to primaries and red to secondaries. Solid lines show the evolu¬ 
tion of rotation including the tidal interactions between the stars. 
The dashed lines show the respective evolution if the stars were 
isolated. Independent age and rotation measurements (see Table 
1) are shown with rotation from photometry (blue error bar) and 
spectroscopy (cyan error bar). The models developed here are in 
reasonable agreement with observed rotational periods (see text for 
discussion). 
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Figure 7. SW and XUV flux evolution for Kepler-16. While harsh XUV conditions exist in the BHZ (top two panels), very favorable SW 
conditions are estimated (bottom two panels) to be better than in the solar system. A Mars-sized planet at the inner edge of the Kepler-16 
BHZ resides in a plasma environment slightly better than Mars. 


The focus of this paper is on the question of how the 
evolution of rotation and activity translates into a vari¬ 
able plasma and radiation environment within the BHZ 
of three well-known circumbinary planets. Both effects 
impact the evolution of potential planets and moons and 
their habitability. 

Figures 7-8, 9 show the XUV radiation and SW fluxes 
as calculated within the BHZ of the binaries (gray shaded 
region). Instantaneous and integrated fluxes are shown. 
In order to evaluate the level of stellar aggression expe¬ 
rienced around each Kepler binary, we compare these re¬ 
sults with those calculated for the SRM (solid thin lines). 
We have compared fluxes around binaries with those for 
Venus, Earth, and Mars in the solar system. For each 
planet we calculate the range of flux experienced by those 
planets for two extreme solar activity evolution models 
(fast and slow initial rotation). 

8.1. Kepler-16 

The conditions around Kepler-16 are shown in Figure 
7. This binary system, which is composed of K and M 
stars, has an XUV radiation environment harsher than 
that of the solar system. Most of the BHZ has been 
exposed to levels comparable to, but mostly larger than, 
those of the solar system planets. Interestingly, the levels 
of XUV radiation received by Kepler-16b and any unob¬ 
served exomoon are almost identical to that received by 
the Earth in the solar system. 


On the other hand, particle fluxes in the BHZ of 
Kepler-16 are almost one order of magnitude lower than 
that solar system HZ levels. We identify the source of this 
significant difference to be the wind acceleration mecha¬ 
nism’s dependence on stellar properties. 

The mass loss rate depends on the amount of heat de¬ 
posited in the photosphere and chromosphere by the tur¬ 
bulent dissipation of Alfvenic waves (Cranmer & Saar 
2011). This heat deposition depends, among many other 
factors, on the third power of the wave amplitude veloc¬ 
ity that goes as with p the photospheric density 

(see Eq. 14 in Cranmer & Saar 2011). In smaller stars, 
with larger photospheric densities, Alfvenic waves will 
propagate with lower amplitudes, and as a result, less 
dissipated heat is available to accelerate the wind. In 
the case of Kepler-16, this, among other less important 
factors, is responsible for a difference, from solar values, 
of two orders of magnitude in the energy available for 
wind acceleration. 

The resulting effect is that, against all odds, the BHZ of 
Kepler-16 (and similar binaries) exhibits enhanced con¬ 
ditions for habitable low-mass planets or exomoons (if 
any). We see in the bottom panels of Figure 7 that 
even a Mars-sized planet at the inner edge of the Kepler- 
16 BHZ will enjoy a plasma environment slightly better 
than Mars. Of course, orbital stability will limit, in this 
case, the very existence of such a planet, but it is still 
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interesting to notice this effect. Moreover, if Kepler-16b 
or a similar planet has a Mars-sized exomoon, then the 
SW appears to pose no serious threat to its atmosphere. 

Is this effect a product of the fact that the planet is 
in a binary system? No, actually. Since the primary is 
just a bit less massive than solar and the companion has 
a very low-mass in a relatively wide orbit, stellar rota¬ 
tion periods have been barely affected by the presence 
of a companion. Therefore, if we locate a planet in a 
scaled orbit around the primary, it will will enjoy similar 
advantageous conditions. 

8.2. Kepler-47 

The case of Kepler-47 is also very interesting (see Fig¬ 
ure 8). In this case, both radiation and plasma environ¬ 
ments, at least with respect to the evolution of secondary 
atmospheres (inset panels), are more planet friendly than 
most of the solar system HZ. The reason, however, is dif¬ 
ferent in Kepler-47 than in Kepler-16. 

First, the primary star is slightly more massive than 
the Sun and produces more XUV radition. It loses mass 
at almost double the solar rate during the very early 
phases of stellar evolution (lower left panel in Figure 8). 
The saturation phase, i.e. the phase during which Rossby 
number is low and chromospheric activity saturates (see 
plateau in XUV flux in the middle panel of Figure 4), 
ends at an early time (20 Myr as compared to 100 Myr 
in the case of the Sun). As a result, the XUV lumi¬ 
nosity decreases much more rapidly during the first cou¬ 
ple of hundreds of Myr. This provides an advantage to 
any protoplanetary and even secondary atmosphere de¬ 
veloped during that phase of stellar and planetary evo¬ 
lution (which extends almost until the age of 1 Gyr). 
In the long run, any planet in the HZ of Kepler-47 will 
have experienced, by the age of the solar system, a cu¬ 
mulative XUV flux lower than most anywhere within the 
Sun’s HZ. 

Around 20 Myr, tidal interaction starts to significantly 
affect primary rotation (see middle panel of Figure 6). 
The effects do not accumulate, however, until the stars 
reach Gyr. 

It is interesting to notice that Kepler-47 offers only 
enhanced conditions in terms of XUV flux early. If our 
assumptions are correct, at around 50 Myr, an Earth-like 
planet in the inner edge of the BHZ will have accumu¬ 
lated XUV flux at a level lower than Venus. In the SW 
case, integrated fluxes during the first several hundreds 
of Myr are enhanced with respect to the Sun. However, 
if secondary atmospheres are common, Mars-sized ob¬ 
jects placed at any stable position within the BHZ could 
preserve their atmospheres. 

8.3. Kepler-453 

The case of Kepler-453 (Figure 9) is close to that of 
the solar system. In terms of insolation, XUV, and SW 
fluxes, the circumbinary system is dominated by the pri¬ 
mary star. The binary period of 27 days and a binary 
eccentricity of 0.05 (see Table I) produce minimal tidal 
interaction for these stars. Differences observed in the 
integrated XUV and SW fluxes arise exclusively from 
differences in stellar properties. 

The most interesting observation concerning the envi¬ 
ronment of Kepler-453 is that even the inner edge of the 


BHZ could be friendly to Venus-like planets (low inte¬ 
grated XUV and SW fluxes). If the fate of Venus habit¬ 
ability was tightly linked to early solar activity, Kepler- 
453 could theoretically harbor more than one Earth-like 
planet within its BHZ. 

9. SUMMARY AND CONCLUSIONS 

We have developed a comprehensive and largely con¬ 
sistent model to calculate the evolution of stellar ac¬ 
tivity in binaries, aimed at constraining the radiation 
and plasma environment of BHZ planets. We apply the 
model to the Kepler binaries with known planets in the 
BHZ, the KBHZ sample, namely Kepler-I6b, Kepler-47c 
and Kepler-453b. 

We find that Kepler-16 probably has a hospitable 
plasma environment for the retention of atmospheres of 
Mars-sized planets and exomoons in the BHZ. Our model 
predicts that the integrated SW flux anywhere in the 
BHZ of the system is lower than that measured at the 
distance of Mars in the solar system. The levels of high- 
energy radiation are, however, similar to those observed 
in the solar system. 

Kepler-47 is the only system among the three in the 
KBHZ sample where tidal interaction has shaped the re¬ 
cent history of stellar activity. High-energy (XUV) radi¬ 
ation levels are affected the most by the evolution of ro¬ 
tation. The historical levels of XUV radiation predicted 
by our model for this system are lower than those mea¬ 
sured in the solar system. This result hints at an intrigu¬ 
ing possibility that Earth-like atmospheres could survive 
desiccation at distances analogous to that of Venus in the 
solar system. The predicted plasma flux in this system is, 
however, too high for Mars-sized planets to retain their 
primary atmospheres. If, however, other mechanisms are 
responsible for the degassing of a secondary atmosphere, 
the late integrated fluxes seem to be lower than those 
predicted in the solar system. 

Einally, the case of Kepler-453 seems to be rather con¬ 
ventional. The circumbinary plasma and radiation envi¬ 
ronment is quite similar to that found in the solar system. 
However, at late times, the integrated high-energy radia¬ 
tion as measured at the inner edge of the HZ (where the 
actual planet Kepler-453b lies) could be better than that 
experienced by Venus. 

These results are not definitive, but they provide a 
solid starting point to consistently assess the key prob¬ 
lem of the radiation and plasma environment of planets 
in binaries. Rotational evolution models of stars in bina¬ 
ries need to be further improved and, more importantly, 
compared with a larger sample of observed rotations in 
binaries. Models of the relationship between stellar ac¬ 
tivity and fundamental properties and rotation period of 
stars, across a wide range of stellar masses, need to be 
further developed and tested. We cannot rely only on 
age-activity relationships found for single-stars in order 
to assess the evolution of activity of stars in binaries. 
Stellar rotation is modified by mutual tidal interaction. 
Any future improvement of these aspects of the mod¬ 
els can be tested using the computational tool we have 
developed and made avaiable at http://bhmcalc.net. 
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Figure 8. SW and XUV flux evolution for Kepler-47. Reduced XUV flux conditions early (top left panel) result in low levels of integrated 
XUV flux (top right panel). However, high SWs early (bottom left panel) result in a Venus-like integrated wind calculation at late times 
(bottom right panel). 
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Figure 9. SW and XUV flux evolution for Kepler-453. Planets residing near the inner edge of the BHZ could be habitable owing to low 
integrated XUV and SW fluxes (right two panels). 
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